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Thermo elastic distortion resulting from optical absorption by transmissive and reflective optics 
can cause unacceptable changes in optical systems that employ high power beams. In advanced- 
generation laser-interferometric gravitational wave detectors, for example, optical absorption is ex¬ 
pected to result in wavefront distortions that would compromise the sensitivity of the detector; thus 
necessitating the use of adaptive thermal compensation. Unfortunately, these systems have long 
thermal time constants and so predictive feed-forward control systems could be required - but the 
finite-element analysis is computationally expensive. We describe here the use of the Betti-Maxwell 
elastodynamic reciprocity theorem to calculate the response of linear elastic bodies (optics) to heat¬ 
ing that has arbitrary spatial distribution. We demonstrate using a simple example, that it can 
yield accurate results in computational times that are significantly less than those required for 
finite-element analyses. 


I. INTRODUCTION 

Advanced-generation interferometric gravitational 
wave detectors, such as Advanced LIGO [T], Advanced 
Virgo [2] and KAGRA [3] are currently being commis¬ 
sioned. Their sensitivity is expected to surpass that 
achieved by first generation instruments by almost an 
order of magnitude in the high frequency region. To 
achieve this, very high circulating power levels (0.5-1 
MW) will be stored within the Fabry-Perot arm cavities. 
At these power levels, even low levels of optical absorp¬ 
tion can lead to significant thermoelastic distortion of 
optical surfaces and unacceptable levels of wavefront dis¬ 
tortion [4], resulting in reduced circulating power and a 
reduction in the efficiency of the detector signal readout. 
Thermally actuated compensation systems will be thus 
used to ameliorate the wavefront distortion. However, 
the thermal time constants for the absorption-induced 
distortion and the compensation are long, typically 12 
hours, and thus incorporating predictive modeling in the 
control systems may prove essential. 

The response of a linear elastic system to heating is 
described by the theory of thermo-elasticity and its ap¬ 
plications to highly symmetric, idealized systems are de¬ 
scribed in many books (see [5] for example). It has 
also been used to develop analytic expressions for less 
idealized optical systems Hi]. The expressions devel¬ 
oped by Hello and Vinet [6] are relevant to the work 
described here, but apply only to cylindrical isotropic 
mirrors heated by coaxial laser beams. 

More complicated systems, which incorporate asym¬ 
metric heating or anisotropic elasticity, can be investi¬ 
gated using finite-element numerical models that apply 
the equations of thermo-elasticity on a three-dimensional 
spatial mesh. For dynamic systems, the thermoelastic 
equations must be solved at each epoch, requiring compu¬ 
tational times that can run to many days. This approach 
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would be untenable for use in predictive feed-forward ac¬ 
tuation to control systems. In such cases, the solution of 
the scalar problem to determine the temperature profile 
throughout the optic can be solved rapidly; the time con¬ 
suming part is solving the tensor-based elasticity problem 
to convert the thermal profile into an elastic distortion. 

The Betti-Maxwell theorem of elastodynamic reci¬ 
procity [7] provides an alternative approach to using 
finite-element methods (FEM) to solve the tensor part 
of the thermoelastic distortion. It has previously been 
used to investigate the excitation of Rayleigh-Lamb elas¬ 
tic waves in a metal plate due to heating produced by a 
line-focused pulsed laser beam assuming that the heating 
is confined to the surface of the plate and it has infinite 
lateral extent 019 ]. In the context of gravitational wave 
detection, it has been used to compute the interferome¬ 
ter’s response to creep events in the fibers that suspend 
the optics HOj. We extend its use to predict thermoelas¬ 
tic distortion of an optic of finite size with asymmetric 
heating. 

We describe here how elastodynamic reciprocity and 
FEM can be combined to provide accurate predictions of 
thermoelastic surface distortion more quickly than using 
FEM alone. In summary, FEM is used to determine the 
response of the optic to a set of ortho normal tractions, or 
pressures —a computationally expensive calculation that 
is performed once for an optic. Then, using reciprocity, 
the distortion due to the instantaneous temperature pro¬ 
file in the optic is calculated using a sum of scalar volume 
integrals that incorporate these responses. The compu¬ 
tational cost of this step is much less than that of a full 
elastostatic FEM evaluation. Additionally, it is amenable 
to parallelization, which would further reduce the com¬ 
putational time. 

The layout of the rest of the paper is as follows: in 
Section H we introduce the Betti-Maxwell theorem of 
elastodynamics and show how it can be used to deter¬ 
mine the surface distortion by careful choice of a suitable 
auxiliary elastic system. We demonstrate its application 
by calculating the distortion of the end face of a cylin¬ 
drical optic that is heated by a Gaussian heat flux that 
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is (a) coaxial with and (b) laterally displaced from the 
axis. The approach and model are described in Sections 
III and IV. Finally, the resulting surface distortions are 
presented in Section V and compared with the results 
of elastostatic FEM calculations. Computation times for 
these two approaches are compared in Section VI 


II. ELASTODYNAMIC RECIPROCITY AND 
THERMAL DISTORTION 

The Betti-Maxwell reciprocity theorem for elastody- 
namics dlS] specifies the relationship between the dis¬ 
placement 'u(r, t) that results from an applied surface 
traction t(r, t) and internal body force /(r, t) for two elas¬ 
tic states of a linear elastic body: 


Consider now applying time-harmonic tractions with 
amplitude tf{rs) = Xn{^s) ^ = 1 , 2 ,.... It is con¬ 
venient to choose Xni'^s) fo be orthonormal, so that 
I Xn{rs)xm{rs)dS = 6„m. Then, expressing the surface 
displacement amplitude as; 

= (4) 

m 

transforms the left-hand term of Eq. § to 

i^dS = / Xn(^ Y amXm(f)dS = a„. (5) 

Therefore 
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where p is the density, u is acceleration, the super¬ 
scripts 1 and 2 represent the two states, and the Einstein 
summation convention is used. If ti(r, t) = and 

/i(r, t) = fi{f)e^^^ then Ui{f^ t) = Ui{r)e^^^^ and thus Eq. 
0 becomes 


an= f afj{f)sfj{f)dV (6) 

Jv 

That is, if the amplitude of the elastic response of the 
optic, , to each of the tractions Xn{^) is known 

then the amplitude of the distortion of the end face of 
the optic, (r) , due to any thermal stress distribution 
can be calculated using Eqs. 0 and 0. 

We shall use this approach to calculate the surface dis¬ 
tortion due to non-uniform heating of a homogeneous 
isotropic body for which 


- t\{^uj{f))dS 

= (// - fi Wi iX) dV (2) 

We shall use this theorem to determine the surface 
displacement (distortion) due to heating of an optic by, 
for example, partial absorption of an incident laser beam. 
Eor the first state, which we shall refer to as the thermal 
state and label T, we assume that the optic is free and 
thus = 0, and there is a non-zero body force due to the 
heating. Since we are interested in the distortion of the 
end face of the optic, we choose the second state, which 
is often referred to as the auxiliary state and we shall 
label A , to have a traction applied to the end face of 
the optic and assume //^ = 0. Thus, Eq. ([^ becomes 

[ tzif)ul{f)dS 

Js 

= f fI{^uf{r)dV = f afj{f^ed{f^dV (3) 
JV Js 

where sfj (r) is the internal strain produced by the trac¬ 
tion tf{r) , and is the internal stress associated 

q^T 

with the body force: ff = • 


^S(r) = r3^AT(f)^,, (7) 

where E is Young’s modulus, a is the coefficient of 
thermal expansion, v is Poisson’s ratio, AT = T{f) — Tq 
, and To is the ambient temperature. Eq. [^thus becomes 

= rS; ~ To)Tr{e^{f)}dV (8) 

III. IMPLEMENTATION 

To determine the distortion of the end-face using reci¬ 
procity, one must first characterize the response of the 
elastic system, sfj (r) , to a set of orthonormal basis trac¬ 
tions = Xn (^) exp [icjt] : n = 1,....,V using an 

elastostatic EEM [11]. 

Zernike functions would be a tempting choice given 
our cylindrical geometry, particularly as they are orthog¬ 
onal to a uniform traction and thus applying the auxil¬ 
iary tractions should not apply net forces to the system. 
However, as shown in Section IV, they are not well suited 
to describing the surface distortion. 

The orthonormal basis tractions we shall use apply a 
non-zero (instantaneous) force to the optic, leading to ill- 
conditioning of the EEM at very low frequencies. We thus 
used a traction frequency of cj = 1 Hz as the response is 
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independent of frequency for frequencies well below the 
first resonance - see [12] for example. 

In all of our numerical tests, we assume a cylindrical 
fused silica optic with height h = 200 mm , radius R = 
170 mm , E = 731 MPa, z/ = 0.17 and a = 0.55 x 10“^ 
K'^ . A radial cross section of the optic and the meshing 
used for the FEM is shown in Fig. 



Figure 1. A radial cross-section of the cylindrical optic, 
showing the mesh used for the FEM. The mesh consisted of 
32000 nodes, and is finest on the heated top surface of the 
test mass. 

We assume heating of the top face by 1 W of power 
absorbed with a Gaussian-distributed fiux: 

Q{x,y) = -^exp[-2{{x-xof + iy- yof)/w'^] 

TTW^ 

where the beam radius w = 53 mm, and radiative cool¬ 
ing of all surfaces of the optic to surroundings at 293 K. 
A thermal FEM m is used to calculate the tempera¬ 
ture distribution, T(r) , resulting from the heating. The 
displacement amplitude for each basis function, , and 
the total displacement, t) , are then calculated us¬ 

ing Eq. I^and Eq. 

IV. CHOICE OF ORTHONORMAL BASIS 
FUNCTIONS 

Choosing a set of ortho normal functions Xn{'^) that 
can describe the surface distortion without requiring a 
large number of functions, which would necessarily in¬ 
clude high spatial frequencies, is crucial as it reduces both 
the number of auxiliary tractions that must be evaluated 
and the requirement for using a fine mesh in the FEM. 
Thus, we describe the choice of basis functions for on-axis 
and off-axis heating of the optic. 

A. Orthonormal basis for on-axis heating 

(xo = 0, yo = 0) 

Zernike polynomials (see Appendix A) are often used 
to describe cylindrically symmetric optical aberrations. 


as they are orthogonal over a circular disc and can be 
normalized. However, as shown in Figure these poly¬ 
nomials are not well suited to describing the distortion. 



Figure 2. Comparison of the surface distortion calculated us¬ 
ing the elastodynamic FEM, Ufem , the sum of the first six 
Zernike components Uz , and the sum of the first six orthonor- 
malized LG components Ulg 

On-axis surface distortion due to the heating can also 
be described using Laguerre-Gauss (LG) functions: 


/2r^\ 


— exp 

2 

V r^ J 

ro J 


where Lp are Laguerre polynomials of order p : 
{0,1,2...} (see Appendix A), r is the radial coordinate 
and ro is a free parameter. These functions are orthogo¬ 
nal only over the infinite plane however. 

Symmetric orthogonalization [12] is therefore used, as 
outlined in Appendix B, to construct linear combina¬ 
tions, Xn : of LG functions that are orthonormal over 
the end face for a given ro . In this type of orthogo¬ 
nalization, the difference between the new and original 
functions is minimized in the least-squares sense [12]. 

The optimum value of ro was chosen as described in 
Appendix C, giving ro = l.bw . The six lowest-order 
orthonomalized-LG functions are defined in Appendix D. 
A comparison of izpem and the sum of these components 
in Fig. shows that the LG basis is much superior to 
the Zernike basis. 


B. Orthonormal basis for off-axis heating 


The distortion due to off-axis heating can be described 
using the sets of functions listed below: 

(a) Hermite-Gauss (HG) functions: 
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where Hi are the (physicists) Hermite polynomials of 
order i : {0,1,2,...} (see Appendix A). These functions 
are orthogonal over the interval x,y : (—oc,oo) . We 
choose rox = = '^o as the heat flux has a circular 

cross section and we shall use xq^i/o « R , and thus 


Zernike polynomial 

an(nm) 

an,FEM(nm) 

Z02 

42.6 

42.9 

Z04 

-15.2 

-15.0 

Z06 

6.3 

5.9 

Z08 

-2.8 

-3.6 


HG^n{x,y) = exp 

(b) Generalized LG functions: 


Table I. Zernike amplitudes calculated usiug reciprocity, an , 
aud thermoelastostatic FEM, an,FEM , for the axisymmetric 
Gaussiau heat flux. 


LG^ir) = Lp 



1 

sin/0 
cos l(j) 


where (j) is the azimuthal angle, and I : {1,2,3,...} for 
p > 0 . We restricted the azimuthal dependence to / = 1 
due to the symmetry of the expected distortion. 

Orthonormalized HG and generalized-LG functions 
were constructed, and an optimized value of tq = lAw 
was selected as discussed above. 

HG functions up to m + n = 15 (136 functions in total) 
were initially used to describe the distortion due to a 
heating beam that was displaced from the center of the 
optic according to (xo,po) =( 05 l 0 mm), (10 mm, 0) and 
(8.7 mm, 5 mm). 

In each case, the distortion was dominated by the same 
17 components, the functions for which are plotted in 
Appendix E. A comparison of i^fem and the sum of the 
dominant 19 components is shown in Fig. [^ . 


V. SURFACE DISTORTION CALCULATED 
USING RECIPROCITY 

We now show how to use the orthonormal bases de¬ 
scribed above with reciprocity to determine the surface 
distortion. In each case, the equilibrium values 

were calculated for the basis tractions and then combined 
with the temperature distribution T(r) from the thermal 
FEM to yield the amplitudes a^. 


A. On-axis heating: Zernike basis 

While Zernike polynomials are not appropriate for de¬ 
scribing the surface distortion in the example presented 
here, they can be used for a reciprocity-based calculation. 
Table I shows a comparison of the reciprocity Zernike am¬ 
plitudes with those calculated by decomposing the distor¬ 
tion predicted by the thermoelastostatic FEM. 



Radial Length [m] 

Eigure 3. Comparison of i^fem and the sum of the 17 dominant 
orthonormalized-HG components for a heat flux offset of 
10 mm 

Orthonormalized generalized-LG functions up to p = 5 
(16 functions in total) were also generated and used to 
describe the distortion due to a heat flux displaced from 
the center of the optic by 10 mm, but they yielded slightly 
poorer agreement with • In addition, since the lower 
order orthonormalized-HG functions appear similar to 
the TEMoi and TEMio eigenmodes observed in optical 
cavities, we chose to use that basis. 


B. On-axis heating: orthonormalized-LG basis 

5 

The UpEu and = J} ^nXn, and the difference be- 

n=0 

tween the two curves are plotted in Fig. Since we are 
not interested in the average displacement of the optics, 
we have set = i^fem at r = 0. The asymmetry of the 
difference is due to non-ideal cylindrical symmetry in the 
FEM meshing. 

C. OfF-axis heating: orthonormalized-HG basis 

19 

The i/pEM and ^nXn and the difference be- 

n=l 

tween the two curves are plotted in Fig. 

Figures and show that even though < 20 auxil¬ 
iary tractions were used to characterize the optic and 
the FEM was restricted to only 30,000 nodes, 

• Elastodynamic reciprocity predicts within < 
1.5% of UpEM over the majority of the incident laser 
beam 

• Displacing the beam by 20% of its radius does not 
degrade the agreement. 
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Figure 4. a) A plot of Ufem and calculated for the on-axis 
heating using the 6 lowest-order orthonormalized-LG func¬ 
tions. b) A plot of Ufem — . 

Additionally, increasing the number of auxiliary tractions 
further improves the agreement, particularly at large ra¬ 
dius. 


VI. COMPARISON OF COMPUTATIONAL 
TIMES 

We compare here the times required to calculate the 
surface distortion using our hybrid FEM-reciprocity ap¬ 
proach and using a conventional thermo-elastic FEM 
analysis. The times are specific to the example of partial 
absorption of a Gaussian-intensity-profile light beam by 
the surface of an isotropic cylindrical optic. 

In both cases we use 32,000 nodes in the EEM calcu¬ 
lations. We have not yet investigated how many nodes 
or auxiliary tractions are required to achieve a particular 
accuracy for each approach, or how this might affect the 
computational times. 

As discussed earlier, our hybrid EEM-reciprocity ap¬ 
proach consists of two parts, the first of which is done 
only once for an optic: 

1. (a) Calculate the elastic response of the optic to 

each of the orthonormal tractions, and store 
these arrays in memory. Here, this consisted 



- 120 ^^-.-■-^.-.--- 

- 0.15 - 0.1 - 0.05 0 0.05 0.1 0.15 

Radial Length [m] 



Figure 5. a) A plot of Ufem and calculated for the off- 
axis heating using the 17 dominant orthonormalized-HG func¬ 
tions. b) A plot of Ufem — • 


of a 32,000 long 6-element array in which the 
3D coordinates and strains at each node were 
recorded for each traction. This part required 
approximately 1 hour per traction. 

(b) Upload 20 responses into memory in prepara¬ 
tion for part 2 required 20 minutes. 

2. At each epoch of interest 

(a) Calculate the thermal induced stress at each 
node using FEM: 90 seconds 

(b) Evaluate the volume integral for each traction 
component using Eq. (6): 3 seconds per trac¬ 
tion. Thus for a serial calculation with 20 trac¬ 
tions, this step required 60 seconds. 

A conventional thermo-elastic EEM calculation for this 
simple problem required about 13 minutes. 

Thus, once the response of the optic has been deter¬ 
mined and uploaded, the hybrid EEM-reciprocity calcu¬ 
lation is at 5.2 times faster is using a serial calculation, 
and 8.7 times faster if using a parallelized calculation of 
the distortion using reciprocity. 
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VII. CONCLUSION 


We have shown how Betti-Maxwell reciprocity can be 
used in combination with thermal finite-element model¬ 
ing to calculate the thermoelastic distortion of a linear 
elastic system. As an example, we described in detail its 
application to calculating the distortion of the end face of 
an isotropic cylindrical glass optic heated by an off-axis 
Gaussian laser beam. Despite using less than 20 auxil¬ 
iary eigenfunction tractions to characterize the optic, the 
distortion calculated using reciprocity agrees to < 1.5 % 
with that calculated using a full thermoelastic FEM over 
the majority of the incident beam. 

The computational time required for the reciprocity 
approach was a factor of 5-8 less than that for the full 
FEM once the optic had been characterized. The ad¬ 
vantage of this approach will thus be most evident in 


cases where the elastic distortion must be calculated fre¬ 
quently, such as in feed-forward control of systems with 
long thermal time constants for example. Parallelization 
of the reciprocity calculation would also allow further im¬ 
provements to the accuracy by employing additional trac¬ 
tions but with negligible additional computational cost. 

Our reciprocity approach can be applied to systems 
with arbitrarily distributed heat fluxes and asymmetric 
anisotropic elastic bodies. Furthermore, while our ex¬ 
ample assumed a free optic, other boundary conditions 
could easily be incorporated into the analysis with an 
appropriate set of auxiliary eigenfunctions. 
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Appendix A: Polynomials used in this paper 


Zernike Polynomial 

orthogonal form with p = r/R) 

Z02 

- 1) 

Z04 

/f (6p^ - + 1) 

Z06 

/^(20p8 - 30p" + 12p2 - 1) 

Z08 

(70/9® - 140(0® + 90/9^ - 20/9^ + 1) 


Table 11. Zernike amplitudes calculated using reciprocityun , 
and thermoelastostatic FEM, an.FEM , for the axisymmetric 
Gaussian heat flux. 


n 

Ln ( 3 ?) 

Hn{x) 

0 

1 

1 

1 

-V + 1 

2x 

2 

(x^ -4x +2)/2 

4x^-2 

3 

(-x^ +9x^ - 18a: + 6)/6 

Sx^ — 12x 

4 

(x^^ - 16x^ + 72x^ - 96x + 24)/24 

IGx^^ - 48x^ + 12 


Table III. Laguerre and Hermite polynomials used. 


Appendix B: Summary of symmetric 
orthogonalization 

The linearly independent LG and HG functions, de¬ 
noted here by //c(r) , were orthonormalized over the end 
face of the mirror using the following process [12]: 

1. Galculate the matrix of inner products of the func¬ 
tions: Mj^i = JJ fk{r)fi{r)dS where the inte- 

end face 

gration was evaluated for the mesh used to export 
the data from the FEM. In this work, the data was 
exported on a Imm-pitch mesh and cropped to fit 
within the circular endface of the optic. 

2. Determine the eigenvalues pk and eigenvectors ux 
of the inner product matrix such that Mj^iUix = 
PxUkx 

3. The orthonormalized functions x(r) are then given 
by Xn = ^ E Uknfk 

'' k 
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Appendix C: Optimization of ro 


The optimum ro was chosen to minimize the mean 
squared difference, weighted by the amplitude of the inci¬ 
dent laser beam, between i^fem and the sum of the selected 
orthonormalised components using: 


/ 5 \ ^ 

// (ufeM - CinXn) 

I face ^ n=0 ^ 


exp 


end faee 


{x-xo^ + jy-yo)^ 


dS 


If exp 


end faee 


7/)2 


dS 


where 



Cin 



'^FEM 


XndS 


end faee 


and a new orthonormal set of functions Xn was gener¬ 
ated for each value of ro . The integrations were evalu¬ 
ated using a square array of pitch 1 mm within the end 
face. 

The variation of this mean-weighted-squared difference 
with ro for the axisymmetric heating (^o = 0) is plotted 
in Figure 1^ 



Figure 6. Plot of the mean-weighted-squared difference be¬ 
tween Ufem and the sum of the first six orthonormalized-LG 
components as a function of ro . 

The variation of this mean-weighted-squared difference 
with ro for orthonormalized-HG functions and off-axis 
heating ^o = 10 mm is plotted in Fig. 


Appendix D: Orthonormalized-LG functions used in 
this paper 


Figure 7. Plot of the mean-weighted-squared difference be¬ 
tween Ufem and the sum of the first six orthonormalized-HG 
components as a function of ro . 


n 

COn 

Cln 

C2n 

CSn 

C4n 

C^n 

0 

-0.97 

-0.25 

-0.06 

-0.01 

-0.001 

-0.0001 

1 

0.24 

-0.86 

-0.43 

-0.13 

-0.028 

-0.003 

2 

0.062 

-0.42 

0.69 

0.55 

0.20 

0.039 

3 

0.016 

-0.15 

0.52 

-0.53 

-0.62 

-0.24 

4 

0.005 

-0.055 

0.26 

-0.58 

0.40 

0.79 

5 

-0.005 

0.07 

-0.4 

1.23 

-2.24 

2.17 


Appendix E: The 17 dominant orthonormalized-HG 
functions 


The Hermite Gauss functions up to order n -1- m = 15 
are orthogonalized using symmetric orthogonaliazation. 
Of these 136 modes, the 17 modes that make the largest 
contribution to describing the deformed surface were se¬ 
lected. These modes are shown in Fig. 



Xn — ConLGo + CinLGi + C2nLG2 

+ CSnLGs + C4nLG4 + C^riLG^ 


Figure 8. The 17 dominant orthonormalized-HG functions 




























































